library(RSocrata)
library(tidyverse)
library(lubridate)
library(tidygeocoder)
library(sf)
library(tmap)
library(spatstat)
source("Funciones_propias_ppp.R", encoding = "UTF-8")
La accidentalidad es uno de los problemas más comunes en el planeta provocando perdidas que pueden ir desde materiales hasta la vida misma, por lo que es de particular interés tratar de comprender como se comporta este fenómeno para evitar ser victimas de él.
La problemática en cuestión tiene la complicación de ser algo de naturaleza aleatoria, sin embargo esto no quiere decir que no exista alguna forma de controlarla pues la experiencia ha mostrado que las diferentes ciudades tienen lugares con mayor concentración de accidentes que otras ubicaciones dentro de la misma ciudad.
Debido a que se está trabajando la ocurrencia de accidentes en una ciudad de interés, la distribución Poisson y los procesos puntuales Poisson resultan ser una herramienta adecuada para tratar de explicar el fenómeno.
Una vez presentado el fenómeno que se trabajará, se selecciona la ciudad de Barranquilla para realizar el estudio.
Los datos se extrajeron de la página web GOV.CO la cual tiene diferentes bases da datos de Colombia abiertas al público. ( Accidentalidad en Barranquilla). La base de datos contiene una gran cantidad de variables, sin embargo solo se escogen aquellas de interés para el patrón puntual. La estructura de los datos es presentada a continuación:
accidentes <- read.socrata("https://www.datos.gov.co/resource/yb9r-2dsi.json") %>%
mutate(fecha = ymd(fecha_accidente)) %>%
select(12, 6:8) %>%
filter(year(fecha) == 2021)
knitr::kable(head(accidentes, 10),
col.names = c("Fecha", "Gravedad del accidentes",
"Clase de accidentes", "Sitio accidente"))
| Fecha | Gravedad del accidentes | Clase de accidentes | Sitio accidente |
|---|---|---|---|
| 2021-01-01 | Solo daños | Choque | CL 31 CR 22 |
| 2021-01-01 | Con heridos | Choque | CL 80B CR 42 |
| 2021-01-01 | Con heridos | Choque | AV DEL RIO ALETA DEL TIBURON |
| 2021-01-02 | Con muertos | Choque | CL 74 21B 23 |
| 2021-01-02 | Con heridos | Choque | CL 72 CR 50 |
| 2021-01-02 | Solo daños | Choque | CR 18 CL 36B |
| 2021-01-02 | Con muertos | Choque | CL 23 CR 15 |
| 2021-01-02 | Con heridos | Choque | VIA 40 CL 85 |
| 2021-01-02 | Solo daños | Choque | CL 72 FRENTE AL 57 68 |
| 2021-01-03 | Solo daños | Choque | CL 45 CR 4 |
Se hace necesario realizarle unos ajustes a la base de datos
considerada puesto que está no incluye las ubicaciones exactas de los
accidentes ya sea en longitud - latitud o proyección UTM, para dicho
propósito se usa la función geo() del paquete
tidygeocoder para convertir las direcciones de los
accidentes en coordenadas de longitud - latitud.
# NO CORRER ESTE CHUNK
direcciones <- paste(accidentes$sitio_exacto_accidente, ", Barranquilla, Colombia", sep = "")
localizaciones <- geo(address = direcciones, method = "arcgis")
write.csv(x = localizaciones, file = "accidentes.csv", row.names = F)
Luego de realizar obtener las coordenadas en longitud - latitud, se
usan múltiples funciones del paquete sf para obtener las
respectivas ubicaciones en formato UTM obteniéndose lo siguiente:
bqlla <- st_read('bqlla.geojson')%>%
st_transform(crs = 3857)
coordenadas <- read.csv("accidentes.csv")
# Se toman 1as primeras 4700 filas para no tener que regenerar el csv
# de localizaciones
accidentes_sf <- cbind(accidentes, coordenadas[8927:13626, 2:3]) %>%
st_as_sf(coords = c('long', 'lat')) %>%
st_set_crs(value = 4326) %>%
st_transform(crs = 3857) %>%
st_intersection(bqlla) %>%
filter(gravedad_accidente == "Con heridos") %>%
select(-fecha, -gravedad_accidente)
knitr::kable(head(accidentes_sf[, c(1:2, 8)], 10))
| clase_accidente | sitio_exacto_accidente | Shape_Area | geometry |
|---|---|---|---|
| Choque | CL 80B CR 42 | 146434812 | POINT (-8329215 1230965) |
| Choque | AV DEL RIO ALETA DEL TIBURON | 146434812 | POINT (-8323586 1231025) |
| Choque | CL 72 CR 50 | 146434812 | POINT (-8327065 1231853) |
| Choque | VIA 40 CL 85 | 146434812 | POINT (-8327354 1235362) |
| Choque | CL 64 CR 66 | 146434812 | POINT (-8325561 1232219) |
| Choque | CR 38 CL 108 | 146434812 | POINT (-8330597 1230153) |
| Choque | CL 82 CR 59B | 146434812 | POINT (-8327880 1233398) |
| Choque | AV CIUDAD CARIBE CL 125 | 146434812 | POINT (-8330059 1230143) |
| Choque | CL 30 CR 6 | 146434812 | POINT (-8323371 1228407) |
| Choque | CL 110 CR 9G | 146434812 | POINT (-8330826 1227455) |
tmap_mode('view')
tm_shape(bqlla)+
tm_polygons(alpha = 0.3, border.alpha = 0.7)+
tm_shape(shp = accidentes_sf)+
tm_dots(size = 0.01)
Uno de los parámetros críticos al momento de modelar un patrón puntual es la intensidad la cual puede ser constante (homogénea) o variable (inhomogénea), por lo tanto es importante iniciar el análisis con una prueba de homogeneidad de la intensidad.
# Definiendo el patron puntual de los datos
datos_ppp <- ppp(x = st_coordinates(accidentes_sf)[, 1],
y = st_coordinates(accidentes_sf)[, 2],
window = as.owin(W = bqlla))
Para iniciar, se divide el mapa de Barranquilla en 25 cuadrantes, si el patrón fuera homogéneo se esperaría que el número de accidentes en cada uno de las subsecciones de la ciudad contenga aproximadamente la misma cantidad de accidentes.
ncuadrantes(datos_ppp)
## x y
## 17 7 2
qc_datos <- quadratcount(datos_ppp, nx = 7, ny = 2)
plot(qc_datos, main = "Accidentes en 2021")
Se aprecia que existe discrepancia entre el número de accidentes por sectores. Claramente al este de la ciudad el número de accidentes es mayor respecto al oeste lo cual es señal de que el patrón puntual en cuestión es de naturaleza inhomogénea.
Adicional a los conteos del número de accidentes en cada uno de los sectores definidos previamente, se usa la prueba \(\chi^2\) para contrastar:
\[ \begin{cases} H_0: \text{La intensidad es constante} \\ H_1: \text{La intensidad no es constante} \end{cases} \]
A continuación se muestra el resultado obtenido
quadrat.test(qc_datos)
##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 638.74, df = 12, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 13 tiles (irregular windows)
Puesto que se tiene un valor - p demasiado pequeño (orden de \(10^{-16}\)) se rechaza la hipótesis de homogeneidad de intensidad, es decir, la intensidad puede ser modelada mediante una relación funcional \(\lambda(x,y)\) y se hace necesario estimarla con algún método ya sea paramétrico o no paramétrico como se vera más adelante.
Previamente se verifico que el número de accidentes en la ciudad tiene intensidad inhomogénea, en esta sección el propósito es estimar dicha función ya sea con metodología paramétrica y no paramética
A continuación se presentan mapas de probabilidad de accidentalidad al realizar una estimación no paramétrica de la función de intensidad basado en funciones kernel los cuales son de la forma:
\[ \hat{\lambda}(x) = \frac{1}{h^2} \sum_{i = 1}^{n} \frac{\kappa\left(\frac{||x-x_i||}{h}\right)}{q(||x||)} \]
El ancho de banda se seleccionó utilizando la función
bw.scott, la cual está implementada para los problemas que
involucra la estimación espacial de la intensidad de puntos observados
en redes lineales, como es el caso de los accidentes viales. Además, se
analizó el error estándar de las estimaciones de densidad obtenido con
esta función comparado con los resultados obtenidos de otras posibles
funciones disponibles en el paquete spatstat y se encontró
que el menor error estándar promedio fue el hallado con la función
seleccionada.
graph_ppp(datos_ppp, round(bw.scott(datos_ppp, isotropic = T), 2),
"Año 2021")
Otra posible solución para estimar la función de intensidad es mediante modelos log - lineales los cuales son de la forma
\[ log \ \lambda (u) = \theta \cdot S(u) \]
donde \(S(u)\) es una función vectorial de valor real evaluada en alguna ubicación \(u = (x, y)\).
Para ajustar estos modelos en R, se usan la función
ppm() del paquete spatstat la cual recibe como
primer argumento un objeto de la clase ppp y como segundo
una formula.
El ajuste de dichos es modelos es sensible a la escala de medición de las ubicaciones; debido a que las ubicaciones se encuentran en metros (los cuales son números muy grandes) es necesario cambiar la escala de medición por una más “pequeña” para evitar problemas de singularidad matricial por lo que se convierten los metros en kilómetros.
# Rescalando los objetos ppp para ajustar los modelos
datos_ppp_km <- rescale(datos_ppp, 1000)
unitname(datos_ppp_km) <- c("km", "kilometers")
Luego de reescalar las ubicaciones se ajustan algunos modelos log-lineales.
# Ajuste de modelos
# Tendencia lineal en x e y
log_linx <- ppm(datos_ppp_km, ~ x)
log_liny <- ppm(datos_ppp_km, ~ y)
log_lin_xy <- ppm(datos_ppp_km, ~ x + y)
## Tendencia cuadratica
log_lin_x_cuad <- ppm(datos_ppp_km, ~ x + I(x^2))
log_lin_y_cuad <- ppm(datos_ppp_km, ~ y + I(y^2))
log_lin_xy_cuad <- ppm(datos_ppp_km, ~ x + I(x^2) + y + I(y^2))
# Tendencia cuadratica en x pero no en y
log_lin_xcuad_y <- ppm(datos_ppp_km, ~ x + I(x^2) + y)
log_lin_inter <- ppm(datos_ppp_km, ~ x*y)
log_lin_cuad_inter <- ppm(datos_ppp_km, ~ x*y + I(x^2) + I(y^2))
\[ \log\ \lambda(x,y) = \beta_0 + \beta_1x \]
par(mfrow = c(1, 2))
plot(predict(log_linx), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_linx, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_2 y \]
par(mfrow = c(1, 2))
plot(predict(log_liny), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_liny, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y \]
par(mfrow = c(1, 2))
plot(predict(log_lin_xy), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_xy, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_{11} x^2 \]
par(mfrow = c(1, 2))
plot(predict(log_lin_x_cuad), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_x_cuad, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_2 y + \beta_{22} y^2 \]
par(mfrow = c(1, 2))
plot(predict(log_lin_y_cuad), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_y_cuad, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_{11} x^2 + \beta_{22} y^2 \]
par(mfrow = c(1, 2))
plot(predict(log_lin_xy_cuad), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_xy_cuad, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_{11} x^2 + \beta_2 y \]
par(mfrow = c(1, 2))
plot(predict(log_lin_xcuad_y), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_xcuad_y, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_{12} xy \]
par(mfrow = c(1, 2))
plot(predict(log_lin_inter), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_inter, which = "smooth",
main = "Residuales suavizados")
\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_{12} xy+ \beta_{11} x^2 + \beta_{22} y^2 \]
par(mfrow = c(1, 2))
plot(predict(log_lin_cuad_inter), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_cuad_inter, which = "smooth",
main = "Residuales suavizados")
Luego de realizar ajustes paramétricos y no paramétricos para la propiedad de primer orden (intensidad), ahora interesa estimar las propiedades de segundo orden las cuales son utiles para determinar la naturaleza del patrón puntual de interés: completamente aleatorio, inhibitorio o agregado.
# patron puntual no parametrico usado previamente
ppp_np <- density.ppp(datos_ppp, sigma = bw.scott,
isotropic = T)
# Implementacion manual de envolventes
# K inhomogena no parametrica
# Kmat_np <- matrix(0, nrow = 513, ncol = 99)
# radios <- seq(0, 3450, length.out = 513)
#
# for(i in 1:99){
# temp <- rpoispp(ppp_np)
# temp <- Kinhom(temp, r = radios, lambda = ppp_np)
# Kmat_np[, i] <- temp$bord.modif
# }
# Knp <- Kinhom(datos_ppp, lambda = ppp_np, r = radios)
# upperK <- apply(Kmat_np, 1, max)
# lowerK <- apply(Kmat_np, 1, min)
#
# Knp_env_data <- data.frame(upper = upperK,
# lower = lowerK,
# observed = Knp$bord.modif,
# theo = Knp$theo,
# r = radios)
# saveRDS(Knp_env_data, "funciones_segundo_orden/Knp_env_data.Rds")
Knp_env_data <- readRDS("funciones_segundo_orden/Knp_env_data.Rds")
# Gráficando las curvas
# k_aux <- data.frame(r = rep(Knp_env_data$r, 2),
# values = c(Knp_env_data$theo,
# Knp_env_data$observed),
# kind = rep(c("Teórica", "Estimada"),
# each = length(Knp_env_data$r)))
#
# Kgg <- ggplot(k_aux, aes(r, y = values, color = kind)) +
# geom_ribbon(aes(ymin = rep(Knp_env_data$lower, 2),
# ymax = rep(Knp_env_data$upper, 2)),
# color = "white", fill = "grey70") +
# geom_path(aes(linetype = kind)) +
# labs(x = "r [metros]",
# y = "K",
# title = "Envolventes para la función K",
# color = "" ,
# linetype = "") +
# theme_minimal() +
# theme(plot.title = element_text(hjust = 0.5)) +
# scale_color_manual(values = c("black", "red"))
#
# Kggplotly <- plotly::ggplotly(Kgg)
# saveRDS(Kggplotly, "Kggplotly.Rds")
Kggplotly <- readRDS("Kggplotly.Rds")
Kggplotly
# Implementacion manual de envolventes
# g inhomogena no parametrica
# gnp <- pcfinhom(datos_ppp, lambda = ppp_np,
# r = radios)
#
# gmat <- matrix(0, nrow = 513, ncol = 99)
#
# for(i in 1:99){
# temp <- rpoispp(ppp_np)
# temp <- pcfinhom(temp, r = radios, lambda = ppp_np)
# gmat[, i] <- temp$trans
# }
#
# upperg <- apply(gmat, 1, max)
# lowerg <- apply(gmat, 1, min)
#
# gnp_env_data <- data.frame(r = radios,
# observed = gnp$trans,
# lower = lowerg,
# upper = upperg,
# theo = 1)
# saveRDS(gnp_env_data, "funciones_segundo_orden/gnp_env_data.Rds")
gnp_env_data <- readRDS("funciones_segundo_orden/gnp_env_data.Rds")
# Gráficando las curvas
# g_aux <- data.frame(r = rep(gnp_env_data$r, 2),
# values = c(gnp_env_data$theo,
# gnp_env_data$observed),
# kind = rep(c("Teórica", "Estimada"),
# each = length(gnp_env_data$r)))
#
# g_gg <- ggplot(g_aux, aes(r, y = values, color = kind)) +
# geom_ribbon(aes(ymin = rep(gnp_env_data$lower, 2),
# ymax = rep(gnp_env_data$upper, 2)),
# color = "white", fill = "grey70") +
# geom_path(aes(linetype = kind)) +
# labs(x = "r [metros]",
# y = "g",
# title = "Envolventes para la función g",
# color = "" ,
# linetype = "") +
# ylim(c(0, 4.4)) +
# theme_minimal() +
# theme(plot.title = element_text(hjust = 0.5)) +
# scale_color_manual(values = c("black", "red"))
#
# g_plotly <- plotly::ggplotly(g_gg)
# saveRDS(g_plotly, "g_plotly.Rds")
g_plotly <- readRDS("g_plotly.Rds")
g_plotly